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We show that the eigenvalue density of a product X = X1X2 ■ ■ ■ Xm of M independent N x N 
Gaussian random matrices in the limit N —> 00 is rotationally symmetric in the complex plane and 
is given by a simple expression p(z,z) = -^a^^i \z\^ 2+ t^ for \z\ < a, and is zero for \z\ > a. The 
parameter a corresponds to the radius of the circular support and is related to the amplitude of 
the Gaussian fluctuations. This form of the eigenvalue density is highly universal. It is identical 
for products of Gaussian Hermitian, non-Hermitian, real or complex random matrices. It does not 
change even if the matrices in the product are taken from different Gaussian ensembles. We present 
a self-contained derivation of this result using a planar diagrammatic technique. Additionally, 
we conjecture that this distribution also holds for any matrices whose elements are independent, 
centered random variables with a finite variance or even more generally for matrices which fulfill 
CN| ' Pastur-Lindeberg's condition. We provide a numerical evidence supporting this conjecture. 
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I. INTRODUCTION 

Initiated by Wigner more than 50 years ago and developed by Dyson, Mehta and others, Random Matrix Theory 
" (RMT) has been successfully applied to various problems ranging from fundamental physics (for a comprehensive 
review see [l|) to engineering and financial applications 0. One of the reasons of such a wide applicability is the 
universality of many results predicted by RMT. Let us take as an example the problem addressed by Wigner, that is 
how to determine the energy spectrum and level spacing distribution of a many-body quantum system. Due to many 
degrees of freedom and sophisticated nature of interactions one has to turn to a statistical description. However, in 
contrast to statistical mechanics where one fixes the Hamiltonian and averages over possible states of the system, 
. Wigner proposed to treat the very Hamiltonian as a random operator, which in turn can be represented as a large 
random matrix. Relevant properties of such a matrix are determined by symmetries of the problem. The great 
discovery of RMT is that many observables are the same for various statistical ensembles of random matrices. 

To illustrate this, let us cite two classical results of RMT. The eigenvalue density of a real symmetric or complex 
Hermitian N x N matrix, whose entries in the upper/lower triangle are independent, identically distributed random 
variables with a finite variance equal to a 2 /N , converges for N — > 00 to a limiting distribution 

o: 

Hj' PW = t — ~y/ 4a 2 — A 2 , for A e [-2a, 2a], (I) 
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known as Wigner's semicircle distribution, one of the best known results of classical RMT. The class of matrices whose 
spectrum converges to the limit law (JlJ is actually much broader and embraces matrices with entries being independent 
random variables which fulfill Pastur-Lindeberg's condition [3|. This is an example of macroscopic universality of 
random matrices. In this paper we concentrate on macroscopic properties and do not discuss microscopic properties 
£T) • of eigenvalue statistics. 

\ An analogous formula for a non-Hermitian random matrix, which is another example of a macroscopic law, reads 

ON . 

O. P(M) = < (2) 




where z = x + iy is a complex number. The distribution ([2]) is called Girko-Ginibre's distribution. The eigenvalue 
density has a rotational symmetry in the complex plane and is uniform inside the circle of radius a. More generally, 
if a matrix has independent but not identically distributed Hermitian and anti-Hermitian degrees of freedom [J], the 
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limit law @ assumes an elliptic form 

!(i-tV<^ for ^(l+r) 2 + <r»(l—r)a - 1 ' 
(3) 
otherwise, 

where a 2 > is an effective scale parameter and r G [— 1, 1] is a flatness of the ellipse. For r = one recovers the 
circular law ©. For r — > ±1 the support of the distribution ([3]) reduces to a cut [—2a, 2a] on the real (for r — > 1) or 
imaginary (for r — \ —1) axis and the distribution itself reduces to a Wigner law ([T]), as one can see by projecting the 
elliptic distribution (J3|) onto the real (imaginary) axis before taking the limit r — > ±1. 

It might be striking that the derivation of the (apparently simple) functional form of p{z, z) for the Girko-Ginibre 
ensemble is less straightforward than the one for the (more complex) Wigner semicircle law. The reason is that there 
are many powerful methods invented for Hcrmitian random matrices: via orthogonal polynomials or Sclbcrg's integral 
Q , supersymmetric method [f| , diagrammatic expansion 0] , Dyson gas @ and free random variables Q . 

In this paper we would like to present a result for non-Hermitian random matrices which is to a large extent 
universal, similarly to the two classical examples cited above. We shall show that the eigenvalue density px(z,z) of 
a product 

X = X 1 X 2 ---X M , (4) 

of M > 2 independent N x N Gaussian matrices for which (Xi^j) = ■■■ = (XM.ij) — and (IXi^I 2 ) = 
crf/N, . . . , (\Xm,vj\ 2 ) = g'm/N for all i,j, assumes in the limit of N — > oo the following form: 

Px(z,z) = < (5) 
[ for \z\ > a, 

where the effective scale parameter a — <?\(Ti . . . &M- This surprisingly simple formula is the main result of our paper. 
What is even more surprising is that this formula holds for a product of independent but not identically distributed 
Gaussian matrices. This means that the individual matrices X^s in the product may come from different Gaussian 
ensembles (GUE, GOE or various elliptic Gaussian non-Hermitian matrices) and the eigenvalue density will always 
be given by (|5|). In other words, even if Xi, . . . , Xm have oblate eigenvalue spectra, with t\ ^ 0, . . . , tm 0, their 
product will have a rotationally-symmetric one. We shall derivate this result with help of a diagrammatic technique 
appropriately tailored to non-Hermitian random matrices [ToL ITTj and to products of random matrices [l2| . In order 
to make the paper self-contained we will also give an introduction to the diagrammatic methods (for a brief review 
see also [Hj]). 

It is tempting to conjecture that the limit law for the product ([5]) holds also for a wider class of matrices, including 
Wigner matrices whose elements are independent, identically distributed random variables with a finite variance or 
more generally, for matrices which fulfill the Pastur-Lindeberg's condition [3]. We will present a numerical support 
for this conjecture. 

The second objective of this paper is to use §5§ in order to verify an interesting conjecture made in Ref. 
saying that if the eigenvalue density p(x, y) of a non-Hermitian matrix X is rotationally symmetric on the complex 
plane z = x + iy, then the marginal distribution p*(x) = J dyp(x,y) obtained by its projection onto the real axis 
or a projection p*(y) = J dxp(x,y) onto the imaginary axis must be equal to the eigenvalue density of the matrix 
(X + X^/y/E or i(X — XT)/vo, respectively, both being Hermitian matrices. If true, this would allow one to calculate 
p(x,y) from p*(x) via the inverse Abel transform. In particular, if one projects the Girko-Ginibre distribution @ 
onto the real (or imaginary) axis, one indeed obtains the Wigner semicircle law: p*{x) = ^-Vc 2 — x 2 , which is the 
same as the eigenvalue density of the matrix (X + X^)/\/8 (or i(X — X^)/y/8). In 14] it was checked numerically 
that the relation seemed to apply also to more complicated ensembles. Here we shall present a counterexample by 
showing that the projection of the eigenvalue density of a product AB of two Hermitian matrices A and B which is 
rotationally symmetric ([S]) is different from the eigenvalue density of the rescaled anti-commutator [AB + BA) /y/8 
and the commutator i(AB — BA)/\/8, so the conjecture is not true. 

II. GENERALITIES 
A. Eigenvalue density and the measure 



We are interested in the eigenvalue distribution of a random matrix X ((4]) being a product of M independent N x N 
real or complex Gaussian matrices. The eigenvalues {Xi} of X are complex since X may be general be non-Hermitian. 
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The eigenvalue distribution is defined by 

Px(*,*) = /^£<y< 3) (*-A«)y (6) 

where z denotes complex conjugate of z. The averaging (...) = J . . . dfi(Xi, . . . , Xm) is done with a factorized 
probability measure, which in the simplest case of identically distributed matrices takes the form 

AT 

dfx(X u X M ) oc H e-T^^D^, (7) 

where DX^ denotes a flat measure. This formula applies to four generic cases of X^ being (a) complex, (b) complex 
Hermitian, (c) real and (d) real symmetric matrices. The parameter a is defined as a = limAr_ J . 00 2Ndof/N 2 where Ndof 
is the number of real degrees of freedom of the matrix X. For (a) the fiat measure is given by DX^ — FJ- dX^^jdX^ij 
or equivalently by DX^ = fj^ d(ReX Mi y)d(ImX /liii ) and a = 4; for (b) DX^ = dX u Y[ i>:j d(ReX Mj y)d(ImX /liii ), 
a — 2; for (c) DX^ — FT.. dX M> j_y, a = 2 and finally for (d) DX^ = J\%>j dX Mj y, a = 1. For (c) and (d) the Hermitian 
conjugate X^ reduces to the transpose Xj. The proportionality symbol in (J7J means that the measure is displayed 
without a normalization constant which is fixed by the condition J dfx(Xi, . . . , Xm) = 1- 

With this choice of a the variance of individual elements (|X Mi ij| 2 ) = 1/N so that the scaling parameters o\ — 
■ ■ ■ = <jm = 1 and hence a — 1 in Eq. §5§ . This means that the eigenvalue density of individual matrices X^ is given 
by the Girko-Ginibre law @ for (a) and (c) and the Wigner law ([T]) for (b) and (d), in both cases with a = 1. For sake 
of simplicity we stick to this choice in the rest of the paper. The spectrum for arbitrary a%, . . . , <jm can be obtained 
by a trivial rescaling. 

Later on we will also consider a general case of matrices from the elliptic ensemble with the eigenvalue distribution 
d3]). We will also consider a product of non-identically distributed matrices, where X±,... ,Xm belong to different 
elliptic ensembles. 



B. The Green's function 



We shall follow here the standard strategy of calculating the eigenvalue density of a random matrix by first calcu- 
lating the Green's function g(z, z) and then using an exact relation between the eigenvalue density and the Green's 
function. Let us recall this relation. Using the following representation of the two-dimensional delta function 



S (2 \z-X) 



one finds 0, 0J-G3 that 



where 



lim —- — 

e-i-0 7T (\Z 



lim --^z 

e->0 7T OZ 



ldg{z,z) 
Px{z,z) = — — . 

7T OZ 
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(10) 



and Ijv is an N x N identity matrix. As we shall see later, the Green's function can be calculated in the limit N — >• oo 
using a summation method for planar Feynman diagrams. It is convenient to think of g(z, z) as a part of a larger 
object OH: a 2N x 2N matrix G with four N x N blocks [13, El: 



G = 




zIn — X ieljsr 
iel N zl N -Xi 



(11) 



Before we continue let us shortly comment on the notation used in the last formula, since we will also use it in the 
remaining part of the paper. The subscripts zz, zz, zz and zz refer to the position of the N x N blocks in the 
corresponding 2N x 2N matrix. In the shorthand notation the arguments {z, z) of a function defined on the complex 



4 



plane are skipped, so the correct reading of, for instance, G zz is G zz = G zz (z, z). We will also use a convention that 
the normalized trace of an N x N matrix denoted by a capital letter will be denoted by the corresponding small letter, 
for instance g zz = jjTiG zz . 

Now coming back to the problem, by inverting the matrix in the brackets on the right-hand side in the last equation 
we can see that the Green's function g(z,z) is equal to the normalized trace of the upper-left sub-matrix, 

g(z, z) = g zz (z, z) = -^Tr G zz (z, z). (12) 

When one calculates the Green's function (fTO]) or the matrix G (JTTJ) , one has to take the limit N — > oo first, and only 
then allow for e —> 0. This comes from the following reasoning. If e = 0, for finite N the function in the brackets 
(. . .) on the right hand side of (JTUJ) has isolated poles on the complex plane. However, in the limit N —> oo the poles 
coalesce and the function becomes non-holomorphic. One cannot then make an analytic continuation of the function 
from holomorphic to nonholomorphic region, as it is done when calculating G by diagrammatic method which utilizes 
0(1/ z) expansion. A small e > is necessary to make G analytic everywhere. If one naively first took the limit e — » 
and only then the limit N — > oo, the matrix G would become block-diagonal: G zz — ((z — X)^ 1 ), G\ z = ((z — X^)^ 1 ) 
and G zz = G zz — 0. However, we shall see that 

q zz (z 1 z) = lim lim / — Tr- tti— — — ; „ / (13) 

y V ' ; e->0JV->oo\iV (zl N - Xi)(zl N - X) + e 2 l N I V ' 

and g zz (z, z) differ from zero in the non-holomorphic region. In Ref. [lO] it was shown that these quantities are related 
to the statistics of left and right eigenvectors of the non-Hcrmitian random matrix ensemble. 

The quantities g zz — g zz are purely imaginary, and 7 = —g zz g zz is a sort of order parameter for non-holomorphic 
behavior, which is positive in a region of the complex plane where the Green's function is non-holomorphic. The effect 
of pole coalescence and the emergence of a non-holomorphic behavior is very similar to the spontaneous breaking of 
a global symmetry in statistical models. In such systems the symmetry is preserved as long as the system size N is 
finite. It may, however, get spontaneously broken in the limit N —¥ 00. Let us take the Ising model as an example. 
Its Hamiltonian is invariant under a global transformation flipping all spins and hence it has a Z2 symmetry. As 
long as the number of spins is finite, the system is ^-symmetric and the average magnetization, which is an order 
parameter, is equal zero. However in the thermodynamic limit, that is when the system size becomes infinite, the Z2 
symmetry gets spontaneously broken below a critical temperature and the average magnetization is non-zero. If one 
first calculated the average magnetization for a finite system and only then took the limit ./V — > 00, the magnetization 
would be zero in this limit for all temperatures. To avoid the problem one can introduce a tiny external magnetic 
field h which weakly breaks the symmetry for finite-size systems. Now, if one first takes the limit N — > 00 and only 
then h—}0, one will obtain the correct result. In our case, the small parameter e plays an analogous role to h and it 
guaranties that non-holomorphic contributions will be correctly picked up for N — > 00. 



C. Linearization 



Let us have a closer look at the function in the brackets in the definition of the Green's function (jlOp . In our original 
problem the matrix X is a product X = X\ . . . Xm of random matrices so it is a non-linear object from the point 
of view of the degrees of freedom that one has to average over. As a consequence the diagrammatic method would 
become very complicated. One can, however, linearize the problem by a trick used in [73 which relies on substituting 
X by a matrix Y of dimensions MN x MN which is linear in Xk's and has eigenvalues closely related to those of X. 
The matrix Y is constructed from A M 's which are placed in a cyclic positions of a sparse MN x MN matrix, 



Y 









\X M 







Xo 



\ 



Xm-i 

/ 



One can immediately discover a relation between eigenvalues of Y and those of X = X\ . 
M-th power Y which gives a block-diagonal matrix 



(14) 



Xm if one calculates the 



Y M = 



Yn 



\ 



(15) 



Y M J 
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with being cyclic permutations of X^s, = X^X^+t . . .X^+m-i (in the cyclic convention X^ + m = X^, and 
= -X"a/). It is easy to see that all blocks Y^ have the same eigenvalues. Indeed, if A is an eigenvalue of Y^ 
to an eigenvector v^, Y^v^ — Xv^, it is also an eigenvalue of Y^_i to the eigenvector = X^-iv^. One can 

see this by multiplying both sides Y^v^ — Xv^ by -X" M _i, obtaining X^-iY^v^ = XX fi-iv^ which is equivalent to 
Y^-iVp—i = Xv^—i. In other words, the matrix Y M has exactly the same eigenvalues as X and each eigenvalue is 
M-fold degenerated. Eigenvalues of X are thus related to those of Y as Xx — Xy ■ The eigenvalue density px{z,z) 
can be calculated from py(w,w) of Y by changing the variables z = w M : 

Quo duj 1 2 

px(z,z) = M-— p Y (w,w) = —\z\- 2+ -p Y (w(z),w(z)). (16) 

The factor M in front of the Jacobian is related to the fact that the transformation z — w M maps the complex plane 
M times onto itself. The problem is thus reduced to finding the spectral density of Y, which is linear with respect 
to Xi, . . . , Xm- The density py(w, w) can be found from the appropriate Green's function. We will show below that 
Py(w,w) is given by a Girko-Ginibre distribution @, irrespectively of M and of n, r% ... and tm- This is a general 
result. In particular, for M = 2 the matrix Y (fT4)l has an anti-diagonal block structure as chiral Gaussian matrices 
which have been intensively studied in the context of spectral properties of the Dirac operator in QCD In this 
case, the form of the eigenvalue density of Y for circular case (t± = T2 = 0) can be inferred from results presented in 
[20l - [22l ] for complex, quaternion real, and real matrices, respectively. 



III. GREEN'S FUNCTION AND PLANAR DIAGRAMS 



In this section we recall the diagrammatic technique of calculating the Green's function. We begin with Hermitian 
matrices and later generalize the method to non-Hermitian ones and eventually to matrices which additionally have 
a block structure like the matrix Y from the previous section. 

Let us make a general comment before we proceed. The diagrammatic method is based on the observation that 
the Green's function G can be interpreted as a generating function for connected two-point Feynman diagrams. In 
the limit N — > 00 only planar diagrams contribute to G since non-planar ones are suppressed by at least a factor 
0(1 /N) [23|, [24|. In this limit one can write a set of two self-consistent algebraic matrix equations which relate G to 
a generating function, E, for one-line irreducible diagrams. The equations are shown schematically in Fig. [T]and will 
be explained later. They can be solved for G. We want to stress that these equations have exactly the same form 
for Hermitian, complex matrices and for matrices with a block structure. They only differ by an algebraic structure 
reflecting indexing of the matrices G and S. 

We finish with a remark that these equations hold for N — > 00. In the context of the discussion about the order of 
taking the limits in (p~3|) this means that one can safely set e = since the limit N — > 00 has already been taken. 



A. Hermitian matrices 



We will first demonstrate the diagrammatic technique on the example of Hermitian matrices and derive the Wigner 
semicircle law (fTJ). Let us assume that A = A\ A = {A a i,}, a = 1, . . . ,N,b = 1, . . . , JV is drawn from an ensemble 
with a probability measure 

dp{A) oce-% TrA2 DA, (17) 

where DA = dA aa Y[ a >b d(ReA a f,)d(ImA a b). The normalization constant, which is implicit in the above formula, 
is fixed by the condition J dp,{A) = 1. The eigenvalues Xi of the matrix A are real. This makes the situation simpler 
than the one for general non-Hermitian matrices discussed in Sec. II. The eigenvalue density can be expressed as [l[ 



P(A) = (^£<K A - A «))' (18) 



where now the delta function is one-dimensional. Also the Green's function G matrix takes a simpler form, 

G = ((Z- A)- 1 ) =({Z- A)- l dp(A). (19) 
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Here Z — zljv, where z is a complex number. The Green's function g(z) = -^TrG(z) is obtained by the Stieltjes 
transform of the eigenvalue density: 



,,( : » = / dX-^r. (20) 



The last equation yields: 



1 



/9(A) = lmg(X + ie), (21) 



7T 



for e — > 0, as follows from a standard representation of the one-dimensional delta function S(x) — — —Jm(x + ie) 1 . 
The above Green's function can be calculated analytically in the large N limit, expanding (fH?|) in terms of powers of 

z-h 

G{z) = Z- 1 + {Z- l A Z- X A Z- 1 ) + {Z- l A Z~ l A Z~ X A Z~ X A Z~ l ) + ... (22) 

Factors Z~ x are independent of A's and thus can be pulled out of the average brackets. What remains are correlation 
functions of the type (A^^ . . . Ai 2n l i 2n ) which by virtue of the Wick theorem can be expressed as products of 
two-point correlation functions (propagators) 

(A ab A cd ) = —5 ad 5 bc . (23) 

This observation allows one to graphically represent equation (I22p as a sum over Feynman diagrams (see for instance 
[25jn. as shown in Fig. [TJ3. Each propagator is represented as a double arc joining two pairs of matrix indices, while 
Z~ b is drawn as a horizontal line joining indices a and b (Fig. [TJ\). In order to calculate G ab one has to sum up 
contributions of all connected diagrams with two external points a, b. For finite N this is not an easy task because 
there are infinitely many diagrams. The problem enormously simplifies in the limit N — > oo since in this limit only 
planar diagrams contribute to the leading term of 1/N expansion and all non-planar diagrams can be neglected 
23, 24]. It turns out that all planar diagrams can be summed up using an old trick known from field theory which 
reduces the problem to a closed set of equations for G. These equations are known as Dyson-Schwinger equations and 
we will discuss them now. 

First, we introduce a generating function E for one-line irreducible diagrams, that is diagrams which cannot be 
split by cutting a single horizontal line (see Fig. HP). Y> ab generates all one- line irreducible diagrams with vertices a 
and b. The two generating functions are related to each other because any diagram from G can be constructed as a 
sandwich of horizontal lines and one-line irreducible diagrams (Fig. [Tp): 

G = Z- 1 + Z- 1 S Z- 1 + Z- 1 Y^Z- 1 S Z- 1 ZZ- 1 + ... = (Z - S) _1 . (24) 

This matrix equation can be viewed as a definition of E. The introduction of E itself does not help to solve the 
problem. However, one can write down an independent equation for E and G. It follows from the observation that 
any one-line irreducible diagram can be obtained from a diagram from G by adding an arc (a propagator) to it 
(Fig. [Tp). This gives 

E Q 6 = G cd —S cd 6 ab = g5 ab , (25) 

or, in matrix notation, E = gtw- Taking trace of both sides we obtain a = g where a = -^TrE is the normalized trace 
of E. The two equations (IM)l and (l2~5j) form a closed set of equations which can be solved for the Green's function 
g(z). Inserting the last equation to with Z = zl^ we have g (z — g) = 1 and hence g{z) = \{z — \/ z 2 — 4) and 
p(A) = 2^V4 — A 2 , as follows from ([2T]) . 

B. Complex matrices 

Let us now discuss how to calculate the Green's function in case of non-Hermitian Gaussian random matrices with 
complex entries (see for instance 13]). The probability measure is now 

dfi(A) cx e - NTrAA1 ]Jd(ReA u )d(W y ), (26) 
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FIG. 1: (A) Feynman rules. (Z -1 )^ is drawn as a line between a and b and the propagator {A a i,A c d} as a double arc joining a 
with d and b with c, respectively. (B) Graphical representation of Eq. (|22|) . The last three displayed graphs correspond to the 
third term in (|22[) . The contribution of the last diagram can be neglected in the large N limit since it is non-planar and has a 
suppressing factor 1/N 2 . (C) Definition of self- energy E. (D) The first Dyson-Schwinger equation which relates G to E. (E) 
The second Dyson-Schwinger equation. 



which corresponds to a = 4 in Eq. ([7]). The propagators are 

(A ab A cd ) = , {A ab A\ d ) = jjSadSbc, 

(27) 

(At b A cd ) = ^6 ad 6 bc , (Al b Al d )=0. 
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It is convenient to think of A and A< as TV x TV sub-matrices of a 2TV x 27V matrix 

^(t£Ho,°0- < 28 > 

The off-diagonal blocks are equal zero for this particular matrix. We use a convention discussed in Section II: the 
position of an TV x TV sub-matrix is denoted by subscripts z, z. We apply the same notation to other 27V x 27V matrices: 
the Green's function, the self-energy X and the matrix Z, 

G =(r zz r")> E =(r"y z ')> z = ( f z Z 7 Z ~ Z ) ■ (29) 

\ "zz <~Tzz J y Lj zz Lj zz J y Zj zz Zj zz J 

Matrix elements of the block G zz of G will be denoted by G a b, elements of G zz by G ab -, etc. In other words, the 
subscripts z and z serve also as templates for the corresponding barred or unbarred indices. For completeness let us 
rewrite the propagators (I27p using this notation: 

(AabAcd) = , (AabA sd ) = jfS ad 8 b5 , 

(30) 

Now we are ready to write down Dyson-Schwinger equations for complex matrices. The first equation is identical 
to Eq. ((24]) . except that now G, X and Z have dimensions 27V x 27V: 



G zz G zz \ / Z zz — 12i zz Z z 

G zz G zz I \ Z zz — Y> zz Z z 



(31) 



This equation is general, but later we will write it for a specific form of Z relevant for the calculation of the eigenvalue 
density. The second equation, which corresponds to (|25[) . can be derived using the propagators defined in Eq. (|30p . 
It can be done separately in each sector zz, zz, zz and zz: 

^•ad = , S a J = j^5 a ^5bcGbc — & a d9zz, 

(32) 

Sad = JjS a dSb c Gbc = &addzz , ^ a d = 0' 

where g zz — -^TtG zz and g zz = jjTrG zz . In matrix notation the last equation can be written as 

' y" y" ) = ( ° ^ ) ■ (33) 
Xzz ^zz J \ gzz^N J X ' 

One should note that the form of this equation is independent of Z while the form of the first Dyson-Schwinger 
equation pip is independent of the propagator structure. If we insert now 

Z = hm ( Z ] N K ] N ) = ( Z l N ° ) (34) 

to Eq. (|3"TT) . remembering that we are allowed to take e — > since all above equations are derived for large TV and 
hence the limit TV — > oo has been taken, we eventually obtain a matrix equation 



( G zz 


G zz \ 




\ G zz 


G zz J 


-( 



z~1_n — S zz — Yi zz 

~^zz zIn ~ Egg 



(35) 



which together with (l33|) forms a closed set of algebraic equations for G(z,z). 

We will now solve this set of equations and then determine p(z,z) using Eq. ©. We first notice that Eq. (|3"3")l 
reduces to a 2 x 2 matrix equation: 

(36) 

where, as before, small letters denote the normalized traces of the corresponding blocks, for instance a zz = -^TrX 2Z . 
Similarly, equation (f3"5j) reduces to 



( Ozz O zz 1 


1 = 1 




9zz \ 


1 &ZZ Ozz ) 




V 9zz 


° ) 



Izz 9zz \ _ ( Z — &ZZ —&zz 
^zz Qzz J V ®zz % ~ @zz 



(37) 
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which, after eliminating cr's with help of Eq. (|3"6]l . leads to 

( 9zz ^ ) = ( z - g - zS Y = ( ~ z 9zs ) . (38) 

\9zz 9zz J \-9zz z J \z\ 2 ~ g zz g zz \9zz z J 

This equation has two solutions. The first one corresponds to g zz = g zz — which gives g zz = z^ 1 and is equivalent to 
the trivial holomorphic solution and hence must be true for large \z\. The second solution corresponds to \z\ 2 — g zz g zz — 
1. In this case the off-diagonal blocks are different from zero and g zz — z. The two solutions match for \z\ 2 = 1. 
Therefore, the first solution holds outside the unit circle and the second one inside the circle. Using the Gauss law 
(O one finds 

{ h for N ^ !. 

p(z,z) = { (39) 
[ for \z\ > 1, 

which is the celebrated Girko-Ginibre distribution (III. H6|. 

To summarize this part, one can write the closed set of algebraic equations for G and E in the large- N limit using 
diagrammatic relations between the generating function for connected two-point planar diagrams (given by G) and 
the generating function for one-line irreducible two-point planar diagrams (given by the free energy E). One can set 
e = in these equations since they are derived already in the limit ./V — > oo. 



C. Complex matrices with a block structure 



We are now ready to calculate the Green's function qy (w, w) for the matrix Y l|14[) which has blocks being 
independent complex non-Hermitian Gaussian matrices [l2[ . The matrix G will be now a 2NM x 2NM matrix having 
four NM x NM blocks G ww , G W w, G^w and which themselves consists of M 2 blocks of size N x N which we 
shall denote by G^, G M p, Gp„ and G^ p respectively, for instance 



Gqinu — 



G n . . . G l 



M 



G M1 . . . G 



(40) 



MM 



There is an analogous block structure for the matrix S. One should distinguish Greek subscripts from Latin subscripts 
giving the position of the matrix elements within the block. For instance, E M p is an N x N sub-matrix of the 
block E wlB and (E^p) t is an element of this sub-matrix. In this convention the normalized trace of a block is 

f/jp = ^TrE M p = -k Yl a =i (^f)oo' *-* ne can now re P ea t the same procedure which we applied to the matrix having 
a single block and derive exact relations between the generating function G and E in the planar limit. The first 
Dyson- Schwinger equation, 



Gww G ww 

Gww Gww 



wInM — ^ww 



WW 

NM — ^ww 



(41) 



is almost identical as ([33)l . except that the blocks and the identity matrices are now of dimensions NM x NM. 
To write the second equation, we need to know the propagators. Let us first define a 2NM x 2NM matrix y, a 
counterpart of A from Eq. 



-y | ^WW ^WW 

^WW ^WW 



Y 
Ft 



(42) 



where Y is cyclic as defined in Eq. (|T4"]) and Y* is anti-cyclic, 

/0 



x\ o 



X M-2 



(43) 
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Since the block matrices J^+i = are assumed to be independent of each other, the only non-zero propagators 
are 

(^12,06^21,0?) = O^aftJ^cd) = ■■■ = (yMl,abym,cd) = Jj S ad 5 bc, (44) 

or in short 

(3^12^21) = (^23^32) = ■ ■ • = (3Wim) = T, (45) 

where the tensor T has elements T a i, c d = jj8 a b5 c d, with indices corresponding to the those of the matrices on the 
left-hand side. If we now insert these propagators to the second Dyson-Schwinger equation, we obtain 

£ M A = (46) 

and E^p = E^, y = for fi v. The problem is symmetric with respect to permutation of the matrices X^, so 
g x i = ... = g M! \-j = g W w in the whole uw-block and similarly in the uw-block. Thus the last equation can be 
compactly written as 

~^ww = gww^-NMi ^ww = <?uhiiljVM> (47) 

where Inm is now the identity NM x NM matrix for the whole block, g wa = -j^jTrG W w and g^w = jjjjTtGu, w . 
Inserting T, ww = E^ = and (|47D to (|4"TT) we see that each block on the right-hand side of (j4"Tj) is proportional to 
the identity matrix. Thus equation (|4ip reduces to a 2 x 2 matrix equation for the normalized traces which play the 
role of proportionality coefficients at the identity matrices, 



/ 9 WW 


Qww 




\ 9 WW 


9 WW J 





-1 

9 WW 

~Qww 



(48) 



This is identical to (|38[) for a complex matrix with a single block discussed in the previous section. In other words, 
the Green's function and hence also the eigenvalue density of the matrix Y does not depend on the number of blocks 
in Y and is given by the Girko-Ginibre law [TH, [l6[ 




Girko-Ginibre spectrum into Eq. (|16[) we finally obtain 



p Y (w,w)=< (49) 



This result is valid also for other matrices considered in Eq. ©, that is for real non-symmetric and Hermitian complex 
matrices, as long as M > 1. It is so because what matters is the structure of propagators only, which is the same for 
all mentioned ensembles. In particular, for M = 2 one can deduce this formula from considerations of chiral ensembles 
[20l - l22l | . In the next section we shall show how to derive the above result for the product of M elliptic complex and / or 
real matrices with different oblateness parameters t\ 7^ . . . 7^ tm ■ Now we will only observe that by inserting the 

T^\z\~ 2+ ^ for|z|<l, 
Px (z,z)=px(\z\) = { (50) 
for \z\ > 1, 

which completes the derivation of our main result. In Figs. [5] and |3] we show a comparison between the above formula 
and the spectrum of X obtained numerically by diagonalization of finite matrices. The agreement is very good. For 
the spectrum of the product of two Hermitian matrices (GUE) shown in the left panel of Fig. [5] we observe a small 
deviation from rotational symmetry manifesting as an accumulation of eigenvalues along the real axis and a depletion 
of eigenvalues in a narrow strip close to this axis. The number of eigenvalues on the axis grows as yN and the width 
of the strip decreases as 1/yN when N — > 00. This effect is almost identical as the one known for real Girko-Ginibre 
matrices [2^, [22|. If one multiplies three or more GUE matrices the effect disappears. A difference between the 
product of two and the product of more than two GUE matrices is that for two the trace TrXi is real whereas for 
three (or more) it is not. In other words, the constraint of the trace to be real introduces a weak spherical symmetry 
breaking of the eigenvalue spectrum. 
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I 




FIG. 2: Plots of px(z,z) for Xi,X2 being two Hermitian matrices (left), two complex matrices (middle), and for X\ being a 
Hermitian and X2 an elliptic random matrix with (f> — n/3 (right). For each case 100 matrices of size N = 100 were generated. 




1 1 


1 1 1 

*% ~ 

v ~~ 


1 1 


\ 

A - 
*'\ 

, 1 ^ 



1.5 



0.5 




1.5 



FIG. 3: Plots of Mtt\z\ 2 ~^t px(\z\) obtained from simulations for various M and matrix sizes N. The theoretical distribution 
(not shown in the figure) which corresponds to (JSJl is a step function = 1 for < \z\ < 1 and zero otherwise. Left: 

X — X1X2 (M — 2) for N = 100 and Xi,X^ taken from the same ensembles as in Fig. [2] black solid line for Hermitian, 
red dotted line for complex, and blue dashed line for Hermitian elliptic matrices. Middle: M = 2, complex matrices of size 
N = 50, 100, 200,400 (black solid, red dotted, green dashed and blue dotted-dashed lines, respectively). To obtain these plots, 
we averaged spectra of 10000, 1000, 1000 and 500 matrices and constructed histograms of absolute values of their eigenvalues. 
Right: N = 200 and M — 2, 3, 4 (black solid, red dotted and blue dashed lines). For each M, 1000 matrices were generated. 



IV. PRODUCT OF ARBITRARY GAUSSIAN MATRICES (ELLIPTIC ENSEMBLES) 

Let us now consider a general class of non- Hermitian random matrices which include as special cases the well known 
examples of Hermitian (GUE), Girko-Ginibre, and anti- Hermitian ensembles. These "elliptic" ensembles were first 
introduced in Q and can be defined as follows. A complex, elliptic matrix X is obtained as a linear combination 
of two identical, independent Hermitian Gaussian matrices A,B: X — cos{<f))A + ism((f>)B, mixed with an arbitrary 
real mixing parameter </>. Since A and B are independent, the corresponding propagators are (A a f,A c d} = jj5 a dfibc, 
(B a b B c d) = jjS a dSbc, and (A a bB c d) = 0. When one changes variables from A and B to X and X* one finds 

(X ab X cd ) = (Xl h x\ d ) = t ■ ±5 ad 5 bc , (X ab x!j - (xl b X cd ) = ^S ad S bc , (51) 
where r = cos(2<^>). The corresponding integration measure for X reads: 

d/i(X) oc exp \- N Y~~~2 ( Tl XX^ - r^Tr (XX + X^X^) J 1 JJ d(ReX y )d(ImAy). (52) 

For = (r = 1) the matrix X is Hermitian, for (p = tt/2 (t = —1) it is anti- Hermitian while for <f> — 7r /4 (t = 0) it 
is isotropic complex. 
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A. Eigenvalue distribution of a single elliptic random matrix 

One can determine the eigenvalue distribution of X using the same methods as in Sec. Ill B. The only difference is 
that the propagators (X a i>X c d) = {X],X cd ) (|5Tj) do not vanish but are proportional to r. This leads to the following 
modification of the first Dyson-Schwinger equation (|5o) : 



<?zz &zz j _ | T 9zz 9zz 
&zz &zz I \ 9zz T~g zz 



while the second one (1371) stays intact: 



9zz 9zz 
9zz 9zz 



(53) 



(54) 



These equations can be solved for g zz . The solution reads 

9zz = < (55) 




for + J^y. < 1, 



where z = x + iy. The non-holomorphic solution matches the holomorphic one on the ellipse. The eigenvalue density 
is i 

1 dg zz J ^T^y for + (T^7F ^ 



7T dz 



otherwise. 



The parameter r is a measure of flattening of the ellipse on which p(z, z) > 0. For r = the last equation reproduces 
the result for non-Hermitian complex matrices. For r — > 1, the ellipse reduces to a cut on the real axis. In order 
to determine the eigenvalue density in this case one should first project the density for r < 1 onto the real axis: 
p*(x) — J dyp(x,y), and then take the limit r — > 1. One recovers the Wigner semicircle law p*(x) = ^:V4 — x 2 , as 
expected. 



B. Eigenvalue distribution of a product of two or more elliptic random matrices 



We are now interested in the eigenvalue density of the product (j4]) where X^s are drawn from a Gaussian ensemble 
with the measure (f5"2"]) . We shall show that the result is again given by Eq. ([5]) and hence exhibits a large degree of 
universality: it does not depend on r and is exactly the same even if each of the matrices X^ is drawn from a Gaussian 
ensemble with a different flattening parameter r M . We will derive ([5]) for X — X\Xi and then make a comment on 
the generalization to M > 2. 

We will use the linearization and calculate first the eigenvalue density of the matrix Y (fT4^) constructed from X\ and 
X2, having the only non- vanishing propagators given by (|51l) with two parameters t\ and t^. As before, first we have 
to determine the propagator structure for the block matrix y (|42[) and then apply it to derive the Dyson-Schwinger 
equation. The matrix 3^ reads 



y = 



Y 

y f 



/ x x \ 



x 2 

X 



(57) 



V x\ / 

The first non- vanishing propagator comes from the correlations between X^s and A^'s, exactly as in Eq. (|45[) : 

CVia^l) = Otoyra) = TT. (58) 
The next one comes from autocorrelations of A M 's (|5ip which are proportional to r, 



(^12^12) = riT, (3^21) = r 2 T, 



(59) 
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FIG. 4: Example of calculation of <7i2 in Eq. (|61[) . We write the second Dyson-Schwinger equation for £12. The only non- 
vanishing propagator is the one between indices 1,2 and 1,2. Taking the trace of both sides of the equation we arrive at 
C12 = n<?2i. 



and the last one from autocorrelations of XJ.'s 



(60) 



Here T denotes again a tensor with elements T a b c d — jjSadSbc, where a, b are indices of the first matrix and c, d of the 
second one on the right-hand sides of the above equations. All other correlations between the blocks of y vanish. We 
can now write two Dyson-Schwinger equations: 



and 



/ cm 012 Oil 0"l2 

021 022 021 °22 

0"ll 012 011 012 

\ 0"21 022 021 022 






T1521 


922 





T2012 








3n 


522 








T1521 





3ii 


T 23l2 






(61) 



/ 3n 312 3n 9i2 

921 922 921 922 

9ii 912 9ii 9i2 

921 922 921 922 



I w 

V - 



- 0~11 

021 

Oil 

0"21 



w 



-012 

- 0"22 

-0"l2 

"022 



-0n 

-0"21 

- 0'n 

-CT21 



W 



-Ol2 
-O22 
-0"l2 
022 



(62) 



WW: u ww 



In the first equation the off-diagonal blocks are the same as in the previous section fi"6")) . The diagonal blocks a 
now depend on r^'s. As an illustration we show in Fig. 0]a graphical representation of the equation for CT12 = tic^i 
which explains the flip of indices. Let us first look for a holomorphic solution, so assume that off-diagonal blocks of 
g vanish: g ww — g ww — 0. In this case the above equations reduce to 



0"11 0"12 
0"21 022 



T132I 

T2312 



311 312 
321 322 



W - (Tn -012 
— CT 2 1 W — (T 2 2 



(63) 



and the corresponding equations for cr^w and g W w being complex conjugate of those above. This gives 



3ii 
321 



312 
322 



w 
-T2312 



-T1321 
10 



(64) 



which has two solutions: one with <?n — i and the other one with gn — w/y/rpr^. We take the first one because 
it has the correct asymptotic behavior for large w. For this solution we have 322 = — and gi2 — 321 = 0. The 
holomorphic solution has to be sewed with the non-holomorphic one so that at the boundary 512 = 1721 = 0. If we 
assume that these elements vanish also inside the non-holomorphic region (and correspondingly g\2 = g^i — 0), then 
the equation (|5T|) reduces to 



011 0"12 Oil a 12 

021 O22 2 J 0~22 

011 O12 CT11 CT12 

02 1 O22 O21 O22 



52 2 

gn 

322 

gn 



(65) 



with vanishing diagonal blocks. This equation is identical to the equation with n = t 2 = and was discussed in 
the previous section. As we know it gives Girko-Ginibre distribution for the matrix Y and hence we obtain ([5]) for 
X = XiX 2 - 
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One can repeat the whole reasoning for a product of more than two matrices. One finds again that the solution 
l/w valid outside the non-holomorphic region corresponds to vanishing blocks g M „ = — for p ^ v and that it 
can be sewed with the non-holomorphic solution for which the blocks also vanish. This gives a ww — a^w and one 
obtains exactly the same equations as for t\ = ■ ■ ■ = tm = 0. Therefore, for M > 2 the eigenvalue distribution of Y is 
also given by the Girko-Ginibre law. This result is universal: the spectrum of X is given by Eq. (J5J) independently of 
whether we multiply two Hermitian matrices, or Hermitian by generic complex, or Hermitian by anti-Hermitian etc. 
The limiting spectrum is always the same and differs only by finite-size effects. 

One can also extend this result to purely real matrices generated from the ensemble with a measure [i| 

d(J,(X) ~ cxp {-y (TrAA T - tTtXX) j JJdXy . (66) 

ij 

The case r = 1 corresponds to symmetric real matrices, r = — 1 to antisymmetric ones, and r = to isotropic real 
matrices. The diagrammatic equations in the limit TV — ¥ oo are exactly the same as before, because the propagators 
have the same structure. 



V. PROJECTION OF THE SPECTRUM OF A COMMUTATOR OF GUE MATRICES 

In this section we show that the conjecture made in [l4| is not true. Let us consider a matrix X — XiX 2 which 
is a product of two Hermitian GUE matrices Xi,X 2 . According to the formula (©, the eigenvalue density of X is 
px(z, z) = 2 k\z\ f° r \ z \ ^ 1 an d zero otherwise. The projection of this function on the real (or imaginary) axis gives 



1 1 + VI - x 2 

p* (x) = - In — , (67 

IT \X\ 

for — 1 < x < 1. According to [II], this result should be equal to the eigenvalue density p+(x) of {X\X 2 + x\x\)j\f% 
or p~(x) of i(X\X 2 — X\x\)/*ft*>. Up to a scaling factor V8, these spectral densities are equal to the spectra of the 
anticommutator {X\,X2\ or the commutator i\X\,X2\, because X\ = X\,X 2 = X\. Moreover, p-(x) = p+(x) as 
follows from the observation that in the limit N — ¥ 00 all the moments of the commutator and the anticommutator 

are the same: Tr ([X 1 ,X 2 } k ) = Tr ({X 1: X 2 } k ) for all A; = 1,2, 

We calculate now the eigenvalue density p+(x) of the rescaled anticommutator {Xi,X 2 }/V8- We define two 
matrices A = (X% + X 2 )/y/2 and B = (X 1 - X 2 )/y/2 which are also mutually independent Hermitian matrices with 
a factorized probability measure 

dp(A, B) cx e- N / 2T * A2 e- N / 2T * B2 DADB. (68) 

We have {X\, X 2 } = A 2 — B 2 . One can use the technique of free random variables [28] to calculate the eigenvalue 
density of A 2 — B 2 since in the limit N —> 00 the matrices A 2 and B 2 represent free random variables. The 
addition law for a sum of free variables is expressed in terms of an i?-transform or equivalently in terms of a Blue's 
function B(z) which is a functional inverse of the Green's function G(B(z)) = z and takes a simple form B a+ t,(z) = 
B a (z) + Bb(z) — z , where a and b are free random variables. In our case a = A 2 , b = —B 2 . The Green's function G a 
of A 2 is a special case of the Green's function for Wishart distribution, while G& for — B 2 corresponds to a reflected 
Wishart spectrum A — >• — A, and hence 



2 

The Blue functions for both cases read 



G a (z) = 1 - G b (z) = -G a (- Z )= - 1 + ^ T ^ . (69) 



and thus 



B « {z) = l^-ry B ^ = WT7y (70) 



1 1 + z 2 

B a+b (z) = B a + B b - - = — -. (71) 



z z(l — z 2 ) 



15 




FIG. 5: Comparison between p+(x) from Eq. (|74|1 (solid line), p*(x) from Eq. (|67|l (dashed line), and numerical simulations 
(circles) for N = 100 (1000 matrices were generated). 



This equation has to be inverted for G a + b (z) which is the Green's function for the anticommutator: 

l + G a+b (z) 2 
Z G a+b {z){l-G a+b {z)iy 



(72) 



which leads to a cubic equation for G a + b {z). The solution which has the correct behavior G a + b {z) —tl/z for large z 
reads 

- , , _ 1 + 3z 2 + (-1 - I8z 2 + 3V3Vz 2 + llz 4 - z 6 ) 2 / 3 

Gr n _|_f)(z) — = — = ; . [Td) 

3z(-l - 18z 2 + Sv^Vz 2 + llz 4 - z 6 )V3 



Taking into account the scaling factor ^/S we finally arrive at 



2/3 



V8 , ^- , V3 1 + 24x 2 - (1 + 144.T 2 - 6 V^Va; 2 + 88x 4 - 64a; 6 ) 
P+ (x)=-^ImG a+b (xV8 + iO+) = ^ i ' . (74) 

n 6n \x\ (1 + 144x 2 - 6vWz 2 + 88s 4 - 64s 6 ) ' 

This is different from (x) from Eq. (|67p . In Fig. [5] we compare both spectral densities and show also results of 
numerical simulations which perfectly agree with (|74p. This falsifies the conjecture that if the spectrum of a non- 
Hermitian matrix is rotationally symmetric, it can be found by solving the symmetrized or antisymmetrized Hermitian 
problem. 



VI. CONCLUSIONS 



The main result of this paper is that the eigenvalue density of a product of large, centered (with zero mean) 
Gaussian matrices assumes a very universal form (J5j) with a single scaling parameter a representing the radius of a 
circular support in the complex plane and related to the amplitude of fluctuations of matrix entries. The matrices in 
the product do not have to be identical and each of them may belong to a different elliptic ensemble. 

Taking into account the universality of the Wigner's semicircle law or the Girko-Ginibre distribution for matrices 
having their entries drawn from independent distributions, it is tempting to conjecture that our result will also hold 
in this setting. Namely, we suppose that the same asymptotic result holds for products of Wigner matrices having 
independent elements drawn from any centered distribution which fulfills Pastur-Lindeberg's condition [3]. To assess 
the validity of this conjecture we performed numerical simulations, assuming various distributions of elements of the 
matrices. The only requirement was that the variance of the distribution was equal to 1/N. We did not observe any 
deviations from ([5]) for short-tailed distributions. In Fig. [5] we show an example for a uniform distribution with zero 
mean and variance 1/N. 

As far as future projects are concerned, it would be interesting to generalize the discussion to the Gaussian symplectic 
ensemble plj and to study microscopic properties of eigenvalues of the product of various types of Gaussian matrices 
from different invariant ensembles 20 22]. It would also be interesting to analytically derive the formula for the 
eigenvalue distribution of the product of M matrices of finite size N (see Fig. [2] in the middle). For the Girko-Ginibre 
ensemble 29] it is given by p(z) = crfc(\/2(|z| — 1)^/N)/(2tt). We expect a qualitatively similar behavior also for the 
product of matrices. 
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FIG. 6: Plots of numerically obtained px(|^|) for Xi, X2 being two symmetric matrices which entries (upper triangle) are taken 
from uniform distribution [—^/3/N, W3/N], for N = 200 and for 1000 matrices generated. Dashed line shows the theoretical 
distribution in the limit N — ¥ oo. 

The discussion presented in this paper holds for Gaussian matrices for which the first moment has zero mean, 
(TrX M ) = 0. It would be interesting to check how it changes when (TrA^) ^ 0. This could be the first step towards a 
generalization of Voiculescu's ^-transform composition rule [3(| for calculating the eigenvalue density of asymptotically 
large matrices representing free random variables, to the case when their product has complex eigenvalues. 
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